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ABSTRACT 

We present a new algorithm for detecting filamentary structure FilFinder. The al¬ 
gorithm uses the techniques of mathematical morphology for filament identification, 
presenting a complementary approach to current algorithms which use matched filter¬ 
ing or critical manifolds. Unlike other methods, FilFinder identihes filaments over a 
wide dynamic range in brightness. We apply the new algorithm to far infrared imaging 
data of dust emission released by the Herschel Gould Belt Survey team. Our prelim¬ 
inary analysis characterizes both filaments and fainter striations. We find a typical 
filament width of 0.09 pc across the sample, but the brightness varies from cloud 
to cloud. Several regions show a bimodal filament brightness distribution, with the 
bright mode (filaments) being an order of magnitude brighter than the faint mode 
(striations). Using the Rolling Hough Transform, we characterize the orientations of 
the striations in the data, finding preferred directions that agree with magnetic held 
direction where data are available. There is a suggestive but noisy correlation between 
typical hlament brightness and literature values of the star formation rates for clouds 
in the Gould Belt. 

Key words: ISM: structure - stars: formation - submillimetre: ISM - techniques: 
image processing 


1 INTRODUCTION 

The structure of molecular clouds (MCs) is thought to guide 
the star formation process, establishing both the rate and 
location of star formation. Unfortunately, the molecular hy¬ 
drogen that makes up the majority of the mass in MCs has 
negligible emission at the typical 10 K temperatures that 


characterize the bulk of MCs (Blitz 19931. Furthermore, 


high dust extinction precludes H 2 absorption line mapping of 


cloud structure over large areas (Spitzer et al.|1973|. Thus, 


studies of MC structure must rely on other tracers. Histori¬ 
cally, the most successful of these methods have been molec¬ 
ular line tracers, notably the rotational lines of CO and its 
isotopologues, infrared dust emission, and dust extinction in 
the near infrared and optical. 

Much of the early work studying cloud structure empha¬ 
sized data from large area mapping of molecular emission 


lines (e.g.. Bally et al. 19871. Since the lines are velocity- 
resolved, the additional velocity dimension provided essen¬ 
tial context for understanding the dynamics of MCs. Unfor¬ 
tunately, the molecular line tracers are biased by both chem¬ 
istry, excitation and radiative transfer effects. These effects 
combine to limit the dynamic range of a given molecular line 
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tracer to a relatively narrow range in density ( Evans|1999 1. 
The past decade has seen a shift in attention toward struc¬ 
ture traced by dust, enabled by the thermal emission maps 
provided by the Spitzer and Herschel missions. Dust emis¬ 
sion provides what appears to be a less biased tracer of cloud 
structure, though variations in dust temperature, abundance 


and emissivity still require careful treatment (e.g., Schnee 
|et al.|[2010[ ). Complementary work has been completed us¬ 
ing dust extinction in the near infrared, though these efforts 


have largely focused on nearby clouds (e.g., Alves et al. 2001 
[Lombardi Alves|200i| . 

Qualitative analyses of these tracers forwarded a rich 
view of the structure in MCs, including the significance 


of filamentary structure (e.g., Lynds 


Elmegreen||1979| [Bally et al.|[T987| [Fa garone et al.||199^ 


1962 


Schneider & 


Nagahama et al.[[1998[). However, a quantitative analysis of 


these data was governed by the statistical and numerical 
techniques adopted and these approaches shaped the inter¬ 
pretation of MC structure. These analyses usually followed 
one of two approaches: statistical measurements of cloud im¬ 
age properties or the generation of object catalogs. The sta¬ 
tistical image characterization includes methods like Fourier 
or wavelet analysis which are usually anchored in the the¬ 
ory of turbulence ( Stutzki et al.[[T998 Langer et al.[[l9^ |. 
These methods are most useful when the measurement can 
be linked directly back to the underlying physical proper- 
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ties of the gas. Statistical methods have consistently em¬ 
phasized the turbulent nature of MCs and recent work has 
emphasized retrieving the underlying properties of turbu¬ 


lence from the observations (e.g., Lazarian & Pogosyan 2000 


2004). In contrast, the object cataloging approach splits the 


structure of MCs into an ensemble of substructures and ex¬ 
plore the properties of that ensemble. The classic two algo¬ 


rithms in the field are Gaussclumps (Kramer et al. 1998) 


and Clumpfind (Williams et al. 1994). These methods are 


designed around exploring the structure seen in molecular 
spectral line data cubes. The algorithms search for objects 
that have low (three-dimensional) aspect ratios, either by 
construction (Gaussclumps) or implicitly given the under¬ 
lying pixel grid (clumpfind). Similar methods developed for 
two-dimensional data have been applied to continuum maps 
of the ISM ( Berlin fc Arnouts|1996 Rosolowsky et al. 20101. 
The view of the molecular ISM derived from these object 
identification approaches is frequently described as clumpy. 
This clumpy view is drawn from the analysis of wide area 
mapping of MCs restricted to using only CO as a tracer 
and mapped at relatively poor quality compared to the cur¬ 
rent state of the art: for example, compare data in|Williams| 


et al. (1994) to Narayanan et al. (2008). This view has been 


emphasized by algorithmic construction and analysis. These 
algorithms adapt poorly to recovering structures with large 
aspect ratios, particularly when the structure is not straight 
in the plane of the sky. Additionally, using CO isotopologues 
as molecular gas tracers limits the dynamic range to roughly 
an order of magnitude in volume density for each tracer 
([E vans IpM] ). This density sensitivity obscures filamentary 
structure with a wide range of densities along a line of sight. 

While large-area molecular line mapping now provides 
data sets with wide dynamical range, the observations are 
relatively time consuming. Compared to CO observations, 
dust emission mapping is fast and relatively high resolution. 
The recent dust emission data, most notably the Herschel 
mapping of dust emission in the far infrared and submillime¬ 
tre (70 to 500 fim) present wide area maps with good dy¬ 
namic range and complementary observational biases with 
respect to molecular line mapping. However, dust maps inte¬ 
grate over the line of sight, potentially mixing regions with 
different physical conditions. With the kinematic informa¬ 
tion from CO line emission, some of these overlaps can be 
resolved though line-of-sight projections still influence the 
interpretation ( Beaumont et al(]|2013 |. 

Initial analysis of dust emission maps has focused on 
the filamentary structure of the MC, specifically noting fea¬ 


tures in the images with large aspect ratios (Andre et al. 
|2010[ [^13[ ). Several qualitative descriptions of the filamen¬ 
tary networks have been proposed including a pervasive hub- 


filament structure (Myers 20091, converging filamentary pat- 
terns ([Peretto et al.||2012b[), or parallel networks ([Busque^ 


et al.|2013|. Filaments appear to be ubiquitous though their 


qualitative interpretation can be subject to human judge¬ 
ment and the tracers used. While the dust emission data 
highlight the two-dimensional column density distributions, 
kinematic information from dense gas tracers can be used 
to trace the gas motion. Several individual filaments and 
hub-filament structures have been analyzed in detail and 
have even been argued to be the fundamental structures in 
the cluster formation process (e.g., Myers||2009 Kirk et al. 
2013b|). 


While providing essential guidance, these single-object 
studies must be paired with broader population studies to 
place their results in context. We require a robust character¬ 
ization of filamentary structure that can be used to quan¬ 
tify the properties of filaments for an ensemble of clouds 
across the ISM. With such a quantification of the ensem¬ 
ble, the population can be checked for variation with phys¬ 
ical parameters and compared to both individual cases and 
simulations of cloud structure. This need for a systematic 
characterization of filament properties has driven the devel¬ 
opment or adaptation of several new cataloging algorithms 
in the study of ISM structure. These algorithms approach 
the problem from different perspectives. 

The DisPerSE algorithm ( Sousbie|2011 1 was originally 
developed for identifying structures in three-dimensional 
cosmological simulations. DisPerSE leverages computational 
Morse theory to identify the critical surfaces and subse¬ 
quently draw smooth one-dimensional manifolds (filaments) 
through a density field. At the heart of the method lies a 
Delaunay triangulation of the domain which serves as a well 
defined set of points where the gradients of the underly¬ 
ing field can be rigorously evaluated. The domain is parti¬ 
tioned into regions defined by these critical points and the 
filamentary structures are identified at the borders of these 
domains. By comparing the robustness of the definition sub¬ 
ject to the variation of the points that define the field, the 
algorithm is made robust in the presence of sampling noise. 
The algorithm works in arbitrary dimension and can be ap¬ 
plied to both two-dimensional dust emission data and three- 
dimensional spectral line data cubes. 

Following a similar approach, the Hi-GAL survey team 
( Molinari et al.|2010| is exploring the use of image gradients 
and the local Hessian matrix as a basis for identifying fila¬ 
mentary structure ( [Schisano et al.|2014[ ). The eigenvalues of 
the Hessian matrix identify the direction of elongation and 
can find those sets of points on ridges for which the curva¬ 
ture is large and negative in one direction and small in the 
other. These sets of points define regions which are aggre¬ 
gated into filaments for which the properties can be charac¬ 
terized. The Hi-GAL approach then processes the Hessian 
map using image processing approaches similar to those pro¬ 
posed in Section below. Since the Hessian matrix can be 
generalized to arbitrary dimension, this approach may be 
generalizable to spectral line data, but development to date 
is focused on the two-dimensional images. 


Salji et al. (20151 also present a Hessian-based approach 


but instead uses a multiscale analysis to enhance filamentary 
structures with a range of widths relative to the filter size. 

Finally, the getfilaments algorithm ( Men’shchikov|2013 ) 
identifies filaments through their persistence in the image 
data across multiple spatial scales. The algorithm progres¬ 
sively smooths and subtracts off signal at different spatial 
scales. The smoothing kernel is a Gaussian and increases in 
size exponentially up to the maximum scale probed. The de¬ 
composition is similar to the continuous form of the a trous 
wavelet algorithm ( Starck fc Murtagh||2006' l. Filaments are 
identified as significant structures with lengths minimally 
changed by the convolutions but widths equal to the con¬ 
volution beam. The algorithm is then applied to the data 
across several wavebands, and filaments are reconstructed 
based on a common appearance in several images. 

In this paper, we present a separate approach to iden- 
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tifying filaments based on the application of mathematical 
morphology ( Shih||2009 1. Mathematical morphology has de¬ 
veloped a suite of operations for digital image processing 
based on the underlying pixel grid of the image, which is usu¬ 
ally rectilinear, and the discrete values of the image data. In 
particular, these operators act on grey-valued images, with 
integer value from 0 to 2*^ — 1 where k is usually 1 (binary 
masks), 8, 16 or 32. We identify regions of emission and 
find their skeletons to identify filamentary structures. We 
then prune these skeletons to produce a robust catalog of 
filaments. 

We present our algorithm below. In Section[^we outline 
the algorithm and apply the algorithm to the Herschel Gould 
Belt Survey (Section ^ to generate a catalog of filamentary 
structures (Section |4L We discuss our results in Section]^ 


2 METHODS 

Identifying and segmenting filamentary structure is a chal¬ 
lenging computational problem because there is little prior 
information about the structure of the features in an image. 
This section describes our method of segmenting filamentary 
structure using intensity data and how the main properties 
of each structure are calculated. 

This new algorithm builds off a rich suite of existing 
methods developed in the field of mathematical morphology, 
which is used across science for image processing. Our treat¬ 
ment relies on a well developed theoretical framework that 
is laid out in, for example, |Shih| ( [2009[ hereafter S09). Such 
methods are well validated and public domain implementa¬ 
tions are readily availably We find that this method is ro¬ 
bust to changes in algorithmic parameters and the approach 
can identify filamentary structure over several orders of mag¬ 
nitude in intensity. We refer to our algorithm as FilFinder, 
it is implemented in Python and is publicly availably 


2.1 Isolating Filamentary Structure 


The algorithm begins with a five-step preprocessing se¬ 
quence that identifies regions in the image containing fil¬ 
amentary structure. These steps are demonstrated in Fig¬ 
ure which shows the application of these steps to sim¬ 
ulated data. The image is drawn from an Enzo simulation 
I O’Shea et al.|[2004 l of the column density in a stirred tur¬ 
bulent box with initial conditions chosen to mimic interstel¬ 
lar clouds and produce filamentary structure. We use simu¬ 
lated data to demonstrate the algorithm since the results can 
be produced free from observational noise and then subse¬ 
quently degraded to establish the influence of observational 
effects. First, we fit a log-normal distribution to the bright¬ 
ness data, identifying the mean /r and standard deviation a 
of the log-intensity. The image is flattened using an arctan¬ 
gent transform: I' = /oarctan(///o) where the normalization 
7o = exp(/r-|-2cr). Since arctan(a;) ~ x for x < 1, the effect of 
the transform is to leave emission below the threshold min¬ 
imally changed but values significantly above the threshold 
asymptotically approach 7r7o/2. This transform suppresses 


^ For example, we rely on http://sclkit-image.org/ 
^ https://github.com/e-koch/FllFinder 


objects that are significantly brighter than the filamentary 
structure within the image (usually bright cores). We note 
that the final mask does not change significantly for small 
changes in Iq. We discuss this effect in Appendix]^ Sec¬ 
ond, the flattened data are smoothed with a Gaussian that 
has a FHWM of 0.05 pc - half the expected width of a fila¬ 
ments found by Arzoumanian et al. ( 2011| . The smoothing 
minimizes variations within the filaments, which eliminates 
spurious elements later in the process, while maintaining 
the large-scale structure. Third, we create a mask of the 
filamentary structure by applying an adaptive threshold to 
the smoothed image. This is the crucial step in creating the 
filament mask. Adaptive thresholding considers each pixel 
in the image in comparison to a patch of the image around 
each pixel and uses a criterion to keep or discard that pixel. 
The criterion used here is that the intensity of the central 
pixel must be greater than the median of the neighborhood 
around it. The local threshold is what allows for filamentary 
structure to be detected over the entire dynamic range in 
brightness within the the image. We choose a patch size of 


0.2 pc, as twice the expected filament width (Arzoumanian 


et al. 20111. We note that the choice of patch size affects 


the filamentary mask shape, but this patch size is not used 
when calculating the filament widths (see Section \2A\ . We 
combine this mask with a globally thresholded mask, which 
removes regions below the noise level in the image. 


These first three steps are sensitive to the chosen nor¬ 
malization value for the arctan transform (7o) and the patch 
size used for the adaptive threshold. The vast majority of 
the structure is recognized for a range of reasonable val¬ 
ues of both parameters, however, the connectivity between 
detected regions can change drastically. We find that op¬ 
timal results are obtained by setting these parameters to 
the values explained above but we explore the effects of 
these choices in Appendix]^ Since mathematical morphol¬ 
ogy methods are based on the geometry of the image pixels, 
some instability can arise when the physical scales that we 
use to estimate the algorithm parameters map to small pixel 
sizes. In particular, we find that for regions where the adap¬ 
tive threshold patch size is smaller than 40 pixels across, 
spurious fragmentation occurs. In these cases, the algorithm 
interpolates the image onto a larger grid for the first 3 steps. 
The returned mask is then regridded back onto the original 
image size. 


The remaining two steps clean small, spurious features 
from the adaptive threshold mask. Objects are removed that 
have an area smaller than the typical smallest filamentary 
object: 57r(0.1pc)^, where 0.1 pc is a typical filament width 
and 5 is the aspect ratio for an object to be considered a 
filament. This definition was chosen empirically after exam¬ 
ining the output of multiple datasets. The vast majority of 
filaments make up a network, and so they deviate signifi¬ 
cantly from the elliptical description of filaments. We found 
that, for the dataset analyzed in this work, regions with 
a smaller area than the threshold led to failed fits to their 
radial profiles (Section |2.4[ ), consistent with these being spu¬ 
rious features. Finally, a median filter with a size of 0.05 pc 
is applied to the mask, which removes small features from 
the edges of objects. This step limits spurious branches in 
the filament identification steps that follow. 
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4 Koch & Rosolowsky 



Figure 1. Steps in the masking process, a) Simulated data with Gaussian noise shown on a logarithmic intensity scale, b) Flattened 
using an arctan transform, c) Smoothed using a Gaussian filter, d) Adaptive thresholding, e) Remove unphysical, small objects. 


2.2 Tracing Filaments using Skeletons 


The combination of these five pre-processing steps estab¬ 
lish the mask used for subsequent analysis. Each structure 
within the mask is reduced to a skeleton using a Medial Axis 
Transform ( Blum|[l967 |. This method reduces an object to 
skeleton form (i.e., one pixel wide) such that it is minimally 
connected. The Medial Axis Transform calculates the posi¬ 
tion of the skeleton by computing the minimum distance of 
each pixel in the object to a pixel outside of it. A skeleton 
pixel is one which maximizes this distance with respect to 
its neighbouring pixels along one direction (e.g., vertical), 
but does not in the other direction (e.g., horizontal). This 
can be thought of as finding the circle with the maximum 
radius that will fit into the object. As a simplified example, 
the skeleton of a circular object would simply be the pixel 
at the centre. If the object were then stretched in the verti¬ 
cal direction, its skeleton would be a vertical line. Thus the 
skeleton pixels are points that are the centres of the maximal 
circles that fit into the object. 


After the Medial Axis Transform, we process the re¬ 
sulting skeleton into the filamentary structures. The main 
steps in this process are illustrated in Figure Initially, 
the skeleton contains many small branches resulting from 
small deviations in the mask boundary that get amplified 
under the transform (Figure top panel). These typical 
features need to be “pruned” leaving the dominant skele¬ 
ton feature behind for characterization. Unfortunately, most 
pruning methods present in the literature (S09) shorten the 
length of the skeleton, which we regard as a scientifically rel¬ 
evant property of the filaments. Instead, we complete a full 
analysis of the unpruned skeleton and use the length infor¬ 
mation in the skeleton to identify branches for pruning. We 
first label each pixel in the skeleton by the number of adja¬ 
cent neighbours. There are three relevant cases: end points 
(1 neighbour), body points (2 neighbours) or intersections 



Figure 2. Top: An un-pruned skeleton with intersection points 
removed. Numerical labels denote branches and alphabetical 
ones represent intersections. Middle: The skeleton converted to 
a graph. Bottom: The final skeleton after the pruning process. 


(^ 3 neighbours). Intersection points are removed from the 
skeleton, returning a collection of branches. 

We exploit the fact that each branch must be com¬ 
prised only of end and body points to calculate the length 
of a branch. Each pixel is connected to at most two neigh¬ 
bours (which do not touch each other) and the length of 
each branch is unambiguously defined. The skeleton is then 
converted to a weighted graph, where the nodes of the graph 
consist of the end and intersection points, and the edges are 
the branches. Edges are weighted using a sum of weights de¬ 
rived from the branch length and the average intensity of the 
image along the branch. Both quantities are normalized by 
the sum over all branches in the skeleton, ensuring a weight 
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between 0 and 1. A shortest-path algorithm is run on the 
graph, which identifies the shortest possible path between 
each pair of end points. From this set of paths, we identify 
the longest path through the skeleton. The sum of branch 
lengths for this longest shortest-path is defined to be the 
length of the filament. 

Any remaining branches which are not part of the 
longest path are now subject to pruning. A branch is pruned 
provided its deletion does not affect the skeleton’s connec¬ 
tivity, its length is below a given threshold, and its average 
intensity is insignificant relative to the rest of the skeleton. 
The length threshold is set to 3 times the beamwidth of the 
image; using a shorter threshold did not effectively prune 
spurious branches. The intensity threshold is set to 0.1 of 
the total intensity along the entire skeleton. The value of 
the intensity threshold was found to be insensitive to moder¬ 
ate changes, however it ensures the removal of long spurious 
features over background regions. 

We note that FilFinder does not prohibit loops in the 
skeleton structure. While this can be an artifact of the me¬ 
dial axis transform, we find that the majority of the loops 
describe the actual pattern of emission. The shortest-path 
algorithm used to derive the filament length is not able to 
handle loops in the graph. To overcome this, we only in¬ 
clude the path through the loop that has a higher weight 
when computing the length. The other portion of the loop, 
provided it does not alter the overall connectivity, is then el¬ 
igible to be pruned. To eliminate loops from spurious holes 
in the mask, we fill in holes that have areas smaller than the 
beam area. 


2.3 Plane Orientation and Curvature 

A novel method for parameterizing the orientation of linear 
structure known as the Rolling Hough Transform (RHT) 
was introduced by [Clark et al.| ( |2014| ). We provide a basic 
explanation of the algorithm here. The RHT is based on the 
family of Hough transforms, which map {x, y) space into (r, 
9) space using the equation for a straight line, 

r = xcosO + ysmO, (1) 


where r is the radius from the origin and 9 is the angle from 
a defined direction ( jPuda fc Hart|[l972| . The RHT differs 
from this original Hough transform definition by restricting 
to r = 0 and defining the origin to be at the centre of a disc of 
a given diameter. Then, 0 = 0 is defined to be in the positive 
y-direction. The transform is performed by measuring the 
response within the disc while varying 9 between 0 and tt 
at each specified pixel. The value of the transform around 
each pixel is the angle at which the transform has the largest 
response. 

We perform the RHT using a disc with a diameter of 
3 times the beamwidth at each pixel within a skeleton. 
The diameter was determined to be the smallest diame¬ 
ter which eliminated significant pixelization effects. Smaller 
discs can become dominated by the axes and diagonal di¬ 
rections (9 = 0°,45°,90°) instead of the larger structure. 
Discs significantly larger than filaments can wash out the 
directional information. We sum the response of the trans¬ 
form from all pixels to arrive at a distribution of angles 
in the plane for the skeleton. Unlike in Clark et al. (20141 
where the expectation value of the RHT is compared to the 



270° 


180° 


Figure 3. The output of the RHT for the skeleton shown in 
Figure]^ Note 29 is plotted to visualize the continuous boundary. 
The median is shown by the solid line and the quartiles by the 
dashed lines. The curvature { 59 ) of the filament is defined the 
interquartile range, or the half the angle between the two dashed 
lines, where the ‘half’ is introduced because we are plotting 29 . 


starlight polarization angle, we do not have a quantity to 
set the zero of the distribution. Instead, we calculate the 
weighted directional mean and define it as the orientation 
of the filament. 



/ 'FiWi sin 29i \ 

\TjiWiCos29i J ’ 


( 2 ) 


where Wi is the normalized value of the transform at 9i such 
that EiWi = 1. Here 9 is defined on [—7r/2,7r/2). We calcu¬ 
late the confidence intervals on the directional mean using 
Equation 3 from Fisher & Lewis (19831, 


{9) ± sin ^ 



( 3 ) 


where Ua is the 2 -score of the two-tail probability. 
The quantity a = Ei cos [2wi (61; — (6^))] is the esti¬ 
mated weighted second trigonometric moment and iF — 
[(EiUi; sin0;)^ -I- {TiiWi cos is the weighted length of the 
vector. 

These quantities make the assumption that the distri¬ 
bution is unimodal, which does not hold for all filaments. Be¬ 
cause of this, and since the confidence intervals are limited 
to the range of 9, we define the curvature S9 conservatively 
to be the interquartile range (IQR, 2 = 0.67). Filaments 
with a multi-modal response have an ill-defined orientation, 
however their dispersion is accounted for by the measure of 
curvature. We show an example of the distribution, angle of 
orientation, and curvature in Figurej^ Like the polarization 
angle of starlight, the orientation of a filament is ambiguous, 
so 9 is confined to a range of tt radians. 


2.4 Width 

We calculate a filament’s width by building a radial pro¬ 
file with respect to the skeleton. We consider a distance 
transform using the skeleton as the reference point. To avoid 
counting pixels multiple times, we calculate a global distance 
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transform for an entire image using all of the skeletons. A 
given pixel is then included in profile calculations only for 
the skeleton to which it is closest. The profiles are limited to 
a maximum distance of 0.3 pc to ensure profiles only account 
for the profile of the filament. We found that this threshold 
did not cut off any features of the profile in question for all 
regions that we analyzed (Section]^. We found that this 
limit was not adequate in many cases, particularly when fil¬ 
aments are closely packed together. To ensure we only fit 
the profile of the filament in question, radial profiles are fur¬ 
ther cut by removing portions where the profile begins to 
increase. 

To derive the width and the surface brightness of that 
filament, we fit a Gaussian with a constant offset (back¬ 
ground) to the radial profiles. The radial profile and Gaus¬ 
sian fit for a filament are shown in Figure]^ We define the 
filament width to be the full-width-half-maximum (FWHM) 
value deconvolved by the beamwidth, as described by |Arzou^ 
[manian et al.| ( [2011[ ). 

Not all filaments showed a Gaussian profile, so to obtain 
some measure of filament width in these cases, we adopt a 
parallel, non-parametric method. We estimate the peak and 
background intensities in the profile using the 5th and 99th 
percentiles, respectively. We then interpolate the profile and 
calculate the distance where the intensity drops to half of 
the peak intensity relative to the background. This value 
reproduces the Gaussian FWHM value for a Gaussian pro¬ 
file but still provides a representative width measurement 
for non-Gaussian filaments. The error is estimated by cal¬ 
culating the range of distance between the 45% and 55% 
percentiles in the profile. 

The non-parametric method is used only in the cases 
where a Gaussian does not return a good fit. We quantify 
a ‘good fit’ as one that has errors that do not exceed the 
absolute value of the parameters and the reduced value 
is below ~ 10. The exact cutoff value used is largely insen¬ 
sitive to determining the quality of near Gaussian profiles. 
Those which have significant, non-random deviations from 
a Gaussian return a significantly higher reduced x^- This 
method is used as the last resort in an attempt to derive 
some useful information from the radial profile. In the results 
presented (Section |^, only 2% of the derived widths use 
the non-parametric approach. The non-parametric method 
is used primarily in two situations: there is excessive crowd¬ 
ing with other filaments leading to a lack of valid pixels 
for the profile, or a bright, compact feature not associated 
with the filamentary structure caused the profile to increase 
at greater distance, leading to a severely truncated profile. 
This latter effect occurs for short filaments, as the smaller 
number of surrounding pixels has a greater effect when av¬ 
eraging within each bin. 


3 DATA 

We apply our new filament extraction method to the obser¬ 
vational data acquired as part of the Herschel Gould Belt 
Survey (HGBS; Andre et al.||2010 |. This program surveyed 
16 nearby molecular clouds using the FAGS and SPIRE in¬ 
struments ( [Griffin et al.|[2010 |. The survey produced maps 
of emission with wavelengths of 70 ^m, 160/rm, 250 /rm, 
350 /rm and 500 /rm, which primarily traces dust emission 




Figure 4. Left: The radial profile of the filament shown to the 
right and in Figure The profile is the result of averaging the 
intensity by binning pixels with respect to their distance from the 
skeleton. Right: Skeleton of the filament plotted over the image. 


in these nearby clouds. The Herschel data have provided a 
new view into the (dust) column density structure in nearby 
clouds and were arguably the main reason for the commu¬ 


nity’s increasing interest in filamentary structure (Andre 
et al.|20T0 1 . We used the preliminary release data available 


from the HGBS survey webpage The individual data have 
been described in a series of publications as given in Table 
which also presents adopted cloud properties. At present, 
the HGBS team has not released a final reference set of im¬ 
ages, and our preliminary results described below should be 
revised once a reference set of data is released. 

The HGBS observed nearly all of the star forming 
clouds in the solar neighbourhood, but recent work has 
called attention to the California-Auriga molecular cloud. 
Despite being as massive as Orion (~ 10® Mq), this cloud 
was long ignored owing to the minimal star formation taking 
place therein (Lada et al. 20091. Since it provides a coun¬ 
terpoint to the vigorous star formation taking place in the 
nearby Orion region, we include SPIRE observations by|Har-| 


vey et al. (20131 in our analysis. This provides some ability 


to check whether the differences in star formation rate are 
associated with differences in the filamentary network. 

Maps from the Herschel SPIRE instrument have a zero- 
point intensity offset introduced in the mapping procedure. 
While the relative calibration within the maps is excellent, 
the absolute calibration requires correcting this constant off¬ 
set. Because of significant overlap in the bandpasses, we use 
the 350 /rm data from the second public Planck data re- 
lease ([Planck Collaboration|2015a|. Following Bernard et al. 


( 2010[ ), we determine the offset by convolving the SPIRE 
data to match the Planck resolution, correlating the data 
position-by-position, and fitting a line to the correlation. We 
find that the relative calibration is excellent as measured by 
the slope of the line: within 10% of unity, with the differ¬ 
ences arising in part because of different bandpass shapes 
between the Planck HFI 857 GHz band and the SPIRE 350 


^ http://www.herschel.fr/cea/gouldbelt/en/Phocea/Vie. 
des_labos/Ast/ast_visu.php?id_ast=66 
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Table 1. Adopted Cloud Properties. The masses are calculated using the 
assumptions from Section[5.2|and the Planck offsets. 


Name 


Reference Distance 

SSFR 

Planck Offset 




(pc) 

(Mq yr-i pc-2) 

(Mq) 

(MJy sr-l) 

Aquila 


1, 3 

260 


1.80 

30 000 

85.5 

California Centre 2, 3 

450 


0.32 

23 000 

9.0 

California East 

2, 3 

450 


0.32 

12 000 

10.1 

California West 

2, 3 

450 


0.32 

5200 

14.7 

Chamaeleon 

4, 5, 3 

170 


3.88 

1500 

-879.1 

IC-5146 


6, 3 

460 


0.38 

6800 

20.7 

Lupus I 


7, 3 

150 


0.37 

900 

14.4 

Orion-A Centre 

8, 9, 10 

400 


0.48 

12 000 

32.6 

Orion-A South 

8, 9, 10 

400 


0.48 

18 000 

35.2 

Orion-B 


11, 10 

400 


0.084 

44 000 

26.2 

Perseus 


12, 3 

235 


1.32 

5100 

23.7 

Pipe 


13, 11 

145 


0.032 

1300 

31.7 

Polaris 


1, 14 

150 


0.0 

1000 

9.3 

Taurus 


15, 16, 17 140 


0.147 

2200 

21.2 
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band. The zero-point offset is determined from the con¬ 
stant offset for the line and is reported in Table We note 
that, because of the adaptive thresholding employed in the 
filament identihcation, most of the results are practically 
unaffected by this correction. 

Since the flux ratio between the 350 fim band and the 
other SPIRE bands is a function of dust emissivity and 
temperature, the Planck referencing is ambiguous for other 
bands. While the filament identification algorithm works 
well on these other images, we focus this preliminary analy¬ 
sis on the 350 /rm data because of the clear route to absolute 
calibration. 


4 RESULTS AND ANALYSIS 


In this section, we present results determined purely from 
the filament analysis of the 350 /rm emission alone (@- 
We then examine those results in the context of prelimi¬ 
nary physical interpretations (©■ These latter analyses are 
not complete because dust emission is a function not only 
gas column density but also dust temperature, gas-to-dust 
ratio, and dust emissivity. These confounding variations of 
dust properties preclude a direct translation between emis¬ 
sion and column density. Multiband analysis of filamentary 
structures hnds that they show an increase in dust tem¬ 


perature with filament radius (e.g., Konyves et al. 20101, 
and dust emissivity shows evidence for varying with differ¬ 


ent density histories (Schnee & Goodman 2005 Schnee et al. 


2006 2008 20101. Other efforts are in progress to provide 


definitive maps of the physical conditions in these clouds 
and we confine our analysis to a restrictive linear scaling 
between dust emission and gas column density. While this 
indicates our interpretations are potentially inaccurate, the 
main purpose of this work is to present a new methodology 


and indicate novel directions for robust analysis, and a full 
treatment of these effects is beyond the scope of the present 
work. 

We analyze the SPIRE 350 fim data for the 14 regions 
described in Table above. The maps are convolved using a 
Gaussian beam to a reach a common physical resolution set 
by the angular size of the 25” SPIRE beam at the maximum 
distance present in the analysis (460 pc), which projects to 
0.056 pc. We also regrid the data to a common pixel scale 
based on the smallest distance of the clouds analyzed (140 
pc). This removes any pixel dependencies, which could skew 
the derived properties of each region. We then apply the 
Hlament finding algorithm as described in Section An ex¬ 
ample of the results of the analysis for the Ghamaeleon re¬ 
gion are shown in Figure |5(b)[ For reference, we also plot 
the filaments identifieqj using the DisPerSE algorithm de¬ 


tifiecpl usi 


scribed in Sousbie ( 201l[ ). We used a persistence cutoff of 12 
in the figure. While lowering the persistence reveals fainter 
features, we found that not all of the faint features found 
by FilFinder were recovered by DisPerSE without also re¬ 
covering spurious features. Furthermore, the returned net¬ 
work was largely connected into one feature despite using 
the network breakdown option in the package. The getfila- 
ments code of Men’shchikov (20131 is not publicly released 
at present. The DisPerSE and FilFinder algorithms repro¬ 
duce the general shapes of the brightest filamentary struc¬ 
tures. The exact paths traced by the two approaches are not 
necessarily identical, but the basic features are well repro¬ 
duced. In addition, Fi/Emder identifies numerous features in 
the faint regions of emission that are undetected using Dis¬ 
PerSE. These features agree well with the faint structures 
in the image, particularly those that have been attributed 


^ Following the tutorial: http://www2.iap.fr/users/sousbi6/ 
web/htinl/index38d8 .htnil?post/regular-grid-’/,3A-f ilaments 
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Figure 5. Left: A comparison of filaments returned by the DisPerSE algorithm ( |Sousbie|2011| l (red) and those returned by our method 
(blue) in the Chamaeleon region. Right: A zoomed in view denoted by the box in the left. In general, the filaments identified using 
FilFinder agree well with those identified by DisPerSE. However, FilFinder recovers many more faint, filamentary features away from 
the brightest regions in the image. 


to magnetic features in clouds (Li et al. 2013 Clark et al. 


2014). 


4.1 The Filament Population 

In Figure we present the joint distributions of filament 
properties for all the regions described in Table We adopt 
notation for filament length (L), width {W), curvature {56), 
orientation {6), average surface brightness (1). In Appendix 
[Al we also argue that our parameter distributions are un¬ 
likely to be a significantly affected by the resolution of the 
350 /im data. Any small effects should be removed by the 
convolution of all data to a common physical resolution and 
projection to a common pixel scale. Of the total 369 fil¬ 
aments found, the radial profile fits failed for just 11 fil¬ 
aments due to the complications discussed in Section |2.4| 
In each of these 11 cases, the width was too small to be 
deconvolved; these are the only filaments with widths that 
could not be deconvolved. Overall, we find properties simi¬ 
lar to previous analyses. Filament brightnesses range mostly 
between 20 and 70 MJy sr“^, and are typically around 25 
MJy sr“^. These values correspond to the peak brightness 
above the background in the fit. The median deconvolved fil¬ 
ament width (FWHM) is 0.09 pc with some spread around 
this width (0.06 and 0.14 pc at the 15th and 85th percentiles, 


respectively), in agreement with Andre et al. (2013). There 
is no significant covariance with width: bright filaments show 
no systematic difference in width compared to fainter fila¬ 
ments. Filament lengths are typically 2.2 pc long with the 
majority ranging between 1.4 and 4.4 pc. In exceptional 
cases, a small number of filaments reach up to 20 pc. Note 
that filament lengths are path lengths and the projected dis¬ 
tance between the ends is necessarily smaller. Furthermore, 
this length does not take into account projection effects due 


to the unknown inclination. Based on the ranges of lengths 
and widths, the typical filament aspect ratio is between 7 
and 9. 


4.2 Effect of Cores on Filament Properties 

As indicated in Section we derive filament properties 
based on the original images and do not attempt prepro¬ 
cessing to remove cores. In this section, we briefly describe 
how the inclusion of cores affects the properties of the fila¬ 


ments. As noted in several studies of the HGBS (Konyves 
jet al.|[2010| [Arzoumanian et al.||20lT| [Andre et al.[[2013p , 
the width of cores corresponds well to previously derived 
widths of filaments (~0.1 pc). Based on a by-eye analysis of 
our results, we find that the skeletons returned by FilFinder 
correspond well to tracing the centre of the cores along the 
filament. We note, however, that the skeletonization process 
used does not guarantee this, leading to the possibility that 
the radial profile is artificially wider due to the offset of a 
bright feature. Very few of these cases were found upon in¬ 
specting the radial profiles of each region. We expect that 
this has a negligible effect on the global properties derived. 

Filaments with significantly higher surface brightness 
have been found to have many cores along their extent (e.g.. 


Konyves et al. 20101. This large number of cores likely causes 


a significant increase in the peak of the filament’s radial 
profile. We expect that this is the greatest effect of including 
the core population while deriving filament properties. Since 
the number of low surface brightness filaments far exceeds 
those of high surface brightness, this effect will inflate the 
positive tail of the surface brightness distribution (Figure 

§• 

Finally, we note the dependency of both of these poten¬ 
tial issues on the size of the filament. Many filaments have 
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Figure 6. Key properties derived from for all filaments in the HGBS and the California Molecular Cloud ( [Harvey et al.|2013| l. Note that 
the width refers to the deconvolved FWHM. The dotted lines indicate the 50, 85 and 99.5 percentiles of the histograms. The contours 
represent the areas containing the 50, 85 and 99.5 percentiles. Dots represent individual data points that do not fall within the lowest 
contour. All filaments are above the minimum length cutoff (5 times the beamwidth, or 0.275 pc for the common projection to a distance 
of 460 pc), and so it is not shown in the figure. The beam width of 0.056 pc (projected at 460 pc) is indicated by the red line. 


skeletons which cover regions over a large dynamic range 
in surface brightness. Often the portions of lower surface 
brightness extend substantially from the brighter regions. 
In these cases, both of the issues presented above are un¬ 
likely to have a large effect on the radial profile due to the 
large number of pixels at lower surface brightness. Future 
work that includes a joint analysis of compact core objects 
and filaments will mitigate these possible biases. 


4.3 Region-based Analysis 

The preceding section aggregates all the regions together 
into a single analysis, but this combines widely different re¬ 


gions of the ISM. We take advantage of this broad survey 
of the Gould Belt region to compare the distributions of 
properties in different regions. 

In a by-eye inspection, the brightness of the filamen¬ 
tary networks appears to vary significantly cloud-by-cloud 
(and in regions across the face of the cloud). We summa¬ 
rize these variations by plotting the distributions of filament 
brightness above the background level in Figure[7] This vio¬ 
lin plot shows amplitude of the probability density functions 
(PDFs) as the full width of the ‘violin.’ The probability den¬ 
sity function is generated from a Gaussian kernel density es¬ 
timate of the observed data, which smooths the PDF. The 
quartiles of the data are shown as lines drawn in the distri- 
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Figure 7. A violin plot showing the background subtracted surface brightness along the filamentary structure. The ‘violin’ shows 
an estimate of the PDF for the surface brightness. These regions show significant multimodal structure in their brightness PDFs and 
significant variations between regions. 


butions. This visualization is similar to stacked histograms 
of the data but highlights the differences in shape between 
the regions. Few of the PDFs are Gaussian and many (no¬ 
tably Orion A South, California East and IC-5146) show 
bimodal structure with the bright mode an order of magni¬ 
tude brighter than the faint mode. Next, the brightness dis¬ 
tributions, while overlapping, are significantly different be¬ 
tween the regions. Aquila is dominated by relatively bright 
filaments whereas the Polaris region is substantially fainter. 
The California molecular cloud shares a common mean and 
width for the distribution across all images analyzed for this 
cloud, but the distribution shape varies over the cloud. Simi¬ 
larly, the two parts of Orion A share a common distribution, 
though the “Centre” of the cloud which contains the main 
ridge of the cloud shows no bright mode, indicating that the 
southern portion of the cloud has more mass in the filamen¬ 
tary structure. The Orion B PDF spans the range seen in 
the Orion A images, but still sits at larger values than the 
lowest mass clouds. 

We estimate the influence of filamentary structure on 
the regions as a whole using two different measurements. 
First, we compare the typical filament brightness to the 
background level estimated in the width. On average, fila¬ 
ment brightnesses are I.SIq'j times that of the broad back¬ 
ground on which they are found. This factor does not show 
significant variation from cloud to cloud across the Gould 
Belt. We also use our model to construct an image of the 
filamentary network to quantify what fraction of cloud emis¬ 


sion in the entire image is found within the filamentary net¬ 
work. We find typical values between 4.5% and 23.9% with 
no clear relation to cloud mass or star formation properties. 
Note that this quantification includes all the emission in the 
image, which may spuriously include foreground and back¬ 
ground emission not associated with the cloud. These values 
reflect only the emission contained within the width of the 
filament above the background level derived in the fits; thus 
emission on scales larger than the filament width is not in¬ 
cluded. These values indicate that the filamentary network, 
while prominent describes a relatively a small fraction of the 
total emission in the cloud. 


We compare the width distributions for the filament 


families using a two-sided Kolomorov-Smirnoff test (Press 
et al.|2002 1. In Figure[^we show the negative log-probability 


that the two distributions are drawn from the same par¬ 
ent distribution. The Figure indicates that there are de¬ 
tectable variations in the filament width distributions be¬ 
tween clouds. The clearest associations between populations 
suggests some variation with either cloud mass or distance, 
since the GMCs are, on average, more distant than the low 
mass clouds. The algorithm shows good robustness with re¬ 
spect to distance (see Appendix E and we carry out all the 
analysis with data convolved to a common physical resolu¬ 
tion. However, there may be some distance effects introduced 
because the more distant objects are physically larger in a 
given image and may sample a wider range of properties. 

Overall, we find good evidence for region-based vari- 
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Figure 8. P-values for two-sided Kolmogorov-Smirnoff tests 
showing the probability of the filament distributions being drawn 
from the same parent distribution for the filament width. The 
greyscale value represents the negative logarithm of the probabil¬ 
ity so substantially different distributions are larger (darker). 


ation among the different nearby star-forming molecnlar 
clouds. Several of the Herschel maps observe the same ob¬ 
ject but have divided the map into different regions based on 
observing requirements. In general, we find little variation 
between images observing different parts of the same object. 


4.4 Filament Orientation 


The Rolling Hough Transform (RHT, ^2.3[ ) provides a novel 
means to assess the orientation of linear features. In our 
case, we use it to search for preferred orientations in the fil¬ 
amentary networks that may be caused by magnetic fields 
or the directions of converging flows. We search for the po¬ 
tentially related signatures of the faint striation populations 
relative to brighter filaments. Here we adopt the nomencla¬ 
ture that a ‘striation’ is a faint filament. Despite the fact 
that some clouds show bimodality between bright and faint 
features, we stress that our analysis in Section [4.1| does not 
show evidence for a universal intensity threshold value that 
separates the filaments into bright and faint populations. As 
such we do not specify a threshold to indicate filaments that 
can be definitively considered striations. We show the ori¬ 
entation distributions for six Gould belt regions in Figure 
[^that exhibit the range of orientation distributions in the 
analyzed clouds. We measure the angular distribution of the 
sets of branches that comprise the filaments rather than the 
distribution for entire skeletons as a whole. While this de¬ 
stroys orientation information for individual filaments, this 
approach better describes the orientation of the region’s en¬ 
tire filamentary structure. 


Chamaeleon - Alves de Oliveira et al. (20141 noted that 


the Chamaeleon Molecular cloud shows two distinct filament 


and striation populations (see Figure 5(a) I. By eye exami¬ 
nation of the images suggests that the populations inter¬ 
sect each other perpendicularly. The RHT results show a 


Aquila 





California East 



Figure 9. Filament orientations for six Gould Belt regions. All 
histograms are normalized to unity and the kernel density estima¬ 
tions (solid black line) are computed over a continuous boundary. 
The East-West direction corresponds to 0°. Some clouds show 
preferred orientations (California West, Pipe, Chamaeleon) and 
others are uniform in the brightness distributions (Aquila). 


significant peak in the direction of the striation population 
which dominates the region by surface area. These striations 


lie in the direction of the large scale magnetic field (Alves 
de Oliveira et aL]|2Q14 |. A small peak can be seen approx¬ 


imately in the direction of the bright filaments. Qualita¬ 
tively, the Chamaeleon region has a bimodal distribution of 
the filamentary structure orientation. We attempted to cap¬ 
ture this behaviour quantitatively through different bright¬ 
ness weighting schemes but were not able to improve on this 
qualitative result. 

Aquila - The Aquila region shows no preferred direction 
for the filamentary structure. This is consistent with a by¬ 
eye evaluation of the data. As a whole, the region shows no 
clear bimodality as we see in Chamaeleon. 

California - Despite being part of the same complex, 
the East and West parts of the California Molecular Cloud 
show some differences in the orientation of the filamentary 
structure. The West component has an aligned striation 
structure, though not clearly perpendicular to the promi¬ 
nent filaments as seen in Chamaeleon. This is reflected in 
the output of the RHT (Figure]^, which shows a prominent 
peak centred near 0 radians. There is no noticeable peak 
corresponding to the bright filaments. The central region of 
California (not shown) reveals a similar RHT output as the 
West component. Conversely, the East component yields a 
more uniform distribution with a wider peak centred near 
— 7 r/ 4 . The uniformity makes it more akin to the Aquila re¬ 
gion than the other regions shown in Figure]^ The striation 
structure in this region may have been disrupted by the con¬ 
centrations of YSOs in the region ( [Harvey et al.|26l3 |. 

Pipe - The Pipe molecular cloud has a unique structure 
possibly attributed to a shock ( Peretto et al.|2012a |. There 
is a clear directional preference centred at 6 = 0 , consistent 
with the shape of the bow-front. The small peak corresponds 
to the direction of the connecting structures located behind 
the bow-front. Unlike other regions where the faint filamen¬ 
tary structure dominates the orientation information, the 
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12 Koch & Rosolowsky 


Pipe region is dominated by the prominent filaments. There 
is a noticeable lack of a striation population that is found in 
most other regions in the Gould Belt. This absence may be 
further evidence of the bow-shock dominating the region. 

Orion-A South - The south region of Orion-A has a pre¬ 
ferred directions at about 9 = —1.4 with less significant peak 
closer to 6^ = —0.4. In comparison with Chamaeleon, these 
peaks are less significant and the distribution is more uni¬ 
form. The peak at —0.4 radians corresponds approximately 
to the direction of the bright hlaments in the region. The 
peak at —1.4 radians is associated with the striations, par¬ 
ticularly near the northern part of the image. 

The orientation histograms presented in Figure are 
not weighted by an intensity value associated with each 
branch. We find that such a weighting scheme causes no 
signihcant change in the histogram’s shape for all regions. 
Since intensity scales with filament mass, this points to a 
clear indication that the faint striation structure as a whole 
contributes a signihcant fraction of the mass. The orienta¬ 
tion of these striations is likely connected to the large-scale 
magnetic held around the clouds. In addition to the work 
of Alves de Oliveira et al. (20141 in Chamaeleon, the align¬ 
ment of the faint linear features in the ISM with the large 
scale magnetic held has also been noted in dust extinction 
( Li et al?]|2013 l and Hi ( Clark et al.||2014 l. The sensitivity 
of the FilFinder algorithm to faint hlamentary structure of¬ 
fers a method to quantitatively identify the signature of the 
magnetic held. 


5 DISCUSSION 

Examination of the hlament populations showed signihcant 
region-to-region variation between the different sections of 
the Gould Belt. Since such clouds also show signihcant vari¬ 
ation in their star formation activity, we explore some em¬ 
pirical correlations between hlaments and star formation. 


5.1 Star Formation and Filamentary Structure 

The scientihc interest in hlaments has largely been driven by 
their close connection to the star formation process. Since 
this work presents a means to identify hlaments across a 
range of molecular clouds, we explored whether the proper¬ 
ties of the hlamentary network were connected to the star 
formation rates (SFRs) of clouds. We have adopted the SFRs 
for these regions described in the literature (see Tablej^ and 
parameterize them in terms of a SFR surface density: Esfr. 
In general, these studies identify the numbers of protostars 
of different classes, and assuming an evolutionary timescale 
for the different classes, infer a SFR. These rates are not 
spatially matched to the SPIRE images we analyze in this 
study. In cases where a single SFR is supplied in the litera¬ 
ture for a region, we apply the rate to all images associated 
with that region. 

In Figure we plot Esfr against the median hlament 
brightness for the different images. This parameterization 
explores connection between the strength of the hlamen¬ 
tary structure in a region and star formation rate. The Fil¬ 
Finder approach explicitly includes a smooth background of 
emission underneath the hlamentary network, so the anal¬ 
ysis specihcally measures the magnitude of the hlamentary 
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Figure 10. The SFR densities (see Table[^ plotted against the 
median of the surface brightness along the filamentary structure. 
There is a suggestive trend of the surface density of star formation 
with typical filament brightness though the Chamaeleon region is 
an obvious outlier. 


substructure rather than the cloud surface density overall. 
There is not a clear correlation given the available data; 
however, with the exception of Chamaeleon, images with 
brighter filamentary structure tend to have higher Esfr. 
The Chamaeleon region is small and may have its star for¬ 
mation rate estimates affected by small number statistics. 
The weak correlation is driven by Aquila, Perseus and Orion 
A all having higher than average filament brightnesses and 
Esfr. 

Since the algorithm also fits a background level, above 
which the filament brightness is measured, we compared 
the filament brightness to the background. This background 
emission could be either the relatively smooth bulk of the 
cloud or true background emission, unassociated with the 
star forming region but projected along the line of sight. We 
find that all clouds have a typical filament brightness that 
is typically 1.3 times larger than the background emission 
on which the filaments are found. The bright filaments can 
have filament-to-background ratios of more than 10. 

We also explored the correlation of other properties of 
the filamentary network with star formation rates but did 
not find any obvious connections. However, these results 
may become clearer with the full release of HGBS data in¬ 
cluding column density and temperature results. If borne out 
by a more thorough analysis, this type of correlation would 
support the evolutionary sequence where filaments form and 
gravitationally fragment ( Andre et al.|2013 l, possibly medi¬ 
ated by magnetic fields (Chen & Ostriker, in preparation). 
In this scenario, the early stages of molecular cloud evolution 
would be parameterized by the amplitude of the filamentary 
network within the cloud. Our method provides a method 
to characterize this network. 


5.2 Filament Stability 

We also explore the distribution of filamentary emission that 
would be subject to classical gravitational instabilities. For 
this initial analysis, we adopt a constant dust emissivity at 
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Figure 11. Average linear mass densities of the filaments for a constant = 15 K. The density is calculated by including all pixels 

within the FWHM of the filament. The distributions in the plot are weighted by the number of pixels in the filaments. The dotted line 
corresponds to the thermal critical mass density of 16 Mq pc“^ 


350 fj,m of 0.073 cm^ g ^ 


from the law of Beckwith et al. 


119901 and the assumption that dependence on frequency 


scales with the index /3 = 2. We adopt a uniform dust tem¬ 
perature of Tdust = 15 K as is characteristic of the more care¬ 


ful analysis by Konyves et al. (20101 and Arzoumanian et al. 


120111. Given these assumptions with gas-to-dust mass ratio 


of 100-to-l, we find Egaa = 1.02 Mq pc "^(/sso/MJy sr ). 


Work from Konyves et al. (20101 shows that the dust tem¬ 


perature decreases toward the centres of filaments, meaning 
that the values derived here for the mass will be lower limits 
on the true surface density. Using the modified blackbody 
model for the spectrum, this difference would be a factor of 
4 for Tdast = 10 K. 

The critical mass per unit lengt h for a nearly isother- 
mal cylinder is Mune.crit = Sc^/G (Inutsuka & Miyama 
1997). Taking Tgaa = 10 K as typical for a molecular cloud 


and the mean particle mass in molecular material to be 
p = 2.4, we find Mune.crit = 16 Mq pc“^. We calculate 
the observed Mune assuming that the filament has a Gaus¬ 
sian profile perpendicular to the filament with a FWHM 
as measured previously: W. We average the measured sur¬ 
face densities over the filament within the FWHM width of 
the filament: (Egas) and compute the linear mass density 
as Miine = (Egaa) W/erf(%/ln 2) = 1.23(Egaa)W. The numer¬ 
ical prefactor arises because of the assumed Gaussian pro¬ 
file and the correction required to account for the wings of 
the Gaussian outside of our integration region. More careful 


treatment of projection effects and direct measurements of 
Cs will be required to assess whether a given filament will be 
subject to gravitational instability. Moreover, the instabil¬ 
ity criterion ignores other forms of support such as turbu¬ 
lence or magnetic field configurations that would make this 
threshold value a lower limit. 


Chen et al. (2015 in preparation) argue that the 


presence of velocity gradients across filaments seen in the 


CARMA Large Area Star formation Survey data (Storm 


et al. 2014) indicates that filament formation occurs in 


a magnetized slab. In such an environment, the magnetic 
critical length is L = Bl{2'KpG^^^). Taking the filaments 
as having an axisymmetric Gaussian density profile gives 
Mine 9 Mq pc-^ Here we have 

taken the typical width of a filament to be 0.1 pc consis¬ 
tent with our analysis and the typical magnetic field to be 
B = 30 fiG as a summary value for the results from the 


Planck Collaboration (2015b I. This value is of order the 
gravitational instability value, but depends on local field 
strength and filament width. These parameters can be re¬ 
lated back to mass densities using relationships between pre- 
and post-shock values ( Chen et al.||2015 I but this relation¬ 
ship requires assumptions about the strength of the shock. 
While notable that the different instability mechanisms pro¬ 
duce similar linear mass densities, the measurements that 
ground the theory of filament formation via magnetohydro¬ 
dynamic processes are difficult to obtain from these data 
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alone. This estimate suggests that the surface densities re¬ 
quired for gravitational collapse will also be able to overcome 
magnetic support as well. 

In Figure m we show a violin plot for the linear mass 
density of filaments in the different clouds analyzed here. 
Overall, the linear mass densities of clouds are below the 
stability criterion, but the massive star forming clouds show 
bright filaments well above the instability threshold. While 
the comparison to the values required for instability is in¬ 
teresting, several factors preclude strong conclusions. Some 
filaments in the other regions may be unstable since lower- 
than-assumed dust temperatures can significantly increase 
the mass densities (by a factor likely between 1 and 4). Fur¬ 
ther, the filaments show significant curvature along their 
lengths and this ‘kinking’ would make the systems more 
unstable. We conclude that many filamentary networks are 
near the threshold of instability to fragmentation. Massive 
molecular clouds such as Aquila, California (Centre) and 
Orion show a few bright filaments that would be unstable 
over many different combinations of assumptions, but the 
lower mass clouds do not have the same significant bright 
filaments. 


5.3 Future Work 


6 SUMMARY 


We have presented FilFinder a new algorithm for measuring 
filaments in astronomical images. FilFinder adopts the tech¬ 
niques of mathematical morphology to identify filamentary 
features in image data. This approach differs from other al¬ 
gorithms, which use wavelets and curvelets {get-filaments) 
or critical manifolds (DisPerSE) or Hessian operators (the 
Hi-GAL survey and Salji et n(]|2015 |. The algorithm relies 
on adaptive thresholding, which creates a mask based on 
local changes in brightness. The Medial Axis Transform re¬ 
duces the signal mask to a skeleton. This skeleton is then 
pruned down to a filamentary network. When compared to 
the filaments identified by DisPerSE, FifFinder identifies the 
same bright filaments but also reliably extracts a population 
of faint filaments (striations) in the peripheries of the molec¬ 
ular clouds. 

Such a morphology-based approach enables new types 
of analyses. We apply this method to study the properties 
of the filamentary networks in molecular clouds found in the 
Gould Belt star forming region. We use the Herschel Gould 


Belt Survey data of Andre et al. (20101 supplemented by 


data on the California molecular cloud from [Harvey et al.| 
(20131. These data provide 350 ^m imaging of these clouds 


at a linear resolution of < 0.056 pc. We compare the filament 
populations between different molecular clouds in the Gould 
Belt. Because of the preliminary nature of the data released 
by the Fderschel Gould Belt Survey team, our results are 
tentative and serve to illustrate the types of analysis enabled 
by the new approach. In this preliminary context, we find: 


This paper presents a first look at filamentary structure 
through the lens of a new algorithm. While some of the 
results show that the new method is promising, future work 
is required to make the outcomes more robust. Using ro¬ 
bust maps of the column density and temperature derived 
from multi-band fitting is the first obvious improvement, 
potentially with links to the Planck data on the largest an¬ 
gular scales. With such maps, we expect an increase in the 
mass densities of the filaments and possibly a decrease in 
width as the centres of the filaments gain more importance 
in the fits. The curvature of the filaments and their lengths 
are unlikely to be significantly altered with more physically 
grounded data. 

A primary advantage of this approach, shared with the 
DisPerSE algorithm, is the ability to extend naturally to 
three dimensions. This extension would allow for filament 
extraction from position-position-velocity (PPV) data cubes 
of spectral line emission. The primary barriers to this inno¬ 
vation lie in the careful treatment of intersection and prun¬ 
ing for three-dimensional skeletons. The algorithm is stabi¬ 
lized by the excellent signal to noise present in the Herschel 
data. Typical millimetre-wave PPV data sets have signifi¬ 
cantly poorer sensitivity but new ALMA data sets will have 
the sensitivity required to stabilize the algorithm. 

Finally, a joint analysis of core and filament populations 
will eliminate some of the concerns that cores are distort¬ 


ing filament properties ((4.21. Furthermore, including back¬ 
ground filament structure in a joint analysis with the cores 
will improve estimates of core size, stability and connection 
with their background environment. With the code publicly 
released, we welcome community contribution to its devel¬ 
opment. 


(i) Filament widths are found here to have a median value 
of W = 0.09 pc across the cloud surveyed. The distribution 
shows a variation in width of a factor of ~ 2. The striations 
extracted here alongside the usual bright filaments have a 
similar distribution of widths. 

(ii) The brightness of a filament is typically 1.3 times 
larger than the local background emission. However, because 
filaments are not pervasive, the network comprises about 
10 % of the intensity of the entire image. 

(iii) The brightness distributions of the filamentary net¬ 
work vary significantly from region to region. The large 
molecular clouds in the sample (Orion A and B, Galifor- 
nia) tend to have brighter filamentary networks than fainter 
clouds. Regions without significant star formation (Polaris) 
show significantly lower median filament brightness. Many 
of the distributions are bimodal with a bright mode an order 
of magnitude larger than the faint mode. The distributions 
of filament curvature and width also show significant varia¬ 
tions across the sample. 

(iv) Adapting the Rolling Hough Transform approach of 
[Clark et al.j ( |2014| ), we have examined the distribution of 
orientations for the filaments. Since the faint filaments dom¬ 
inate this analysis by number, we can identify a preferred 
direction in the striation population for some of the clouds 
(California, Chamaeleon, and the Pipe). There is weak ev¬ 
idence for a perpendicular population from the bright fila¬ 
ments in Chamaeleon, consistent with the by-eye analysis of 


Alves de Oliveira et al. (20141. 


(v) We find a suggestive but noisy correlation between 
the brightness of the filamentary network and the liter¬ 
ature estimates of the star formation rate in a region. 
The Chamaeleon region does not agree with this sugges- 
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tive trend, showing a high star formation rate with respect 
to its filament brightness. 

(vi) We examined the stability of the filaments under a 
linear mapping from intensity to column density. The bright 
filaments in the clouds are close to or above the threshold 
for instability. However, the unmeasured quantities in the 
instability criteria and our simple linear mapping surface 
density make both the data and the threshold uncertain by 
a factor of two. 

While such results will need to be re-examined in the con¬ 
text of final data from the Herschel Gould Belt Survey, 
they offer some suggestive directions for the analysis of fil¬ 
aments in the context of far infrared emission data. The 
algorithm can readily be generalized to higher dimensional 
data, enabling a parallel analysis in high-sensitivity position- 
position-velocity data. The code and documentation have 
been publicly released. 
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APPENDIX A: PARAMETER SENSITIVITY 

The creation of the filament mask (Section]^ is dependent 
on two parameters: the size of the patch used for the adap¬ 
tive threshold, and normalization value (7o) used in the arc- 
tan transform. The choice of the parameters does not greatly 
affect the initial detected regions. However, the connectivity 
between these structure may change drastically. We demon¬ 
strate this in Figure |A1[ where the filamentary structure 
in the image from Figure is shown for different combi¬ 
nations of parameters. Choosing a low normalization value 
causes the algorithm to ignore faint connections. Choosing 
a small patch size leads to fragmented regions, which do not 
adequately describe the large-scale structure. Conversely, a 
large patch size will lead to ignoring substructure within fila¬ 
mentary structures. The patch size influences the final mask 
more than the choice of normalization, in fact the mask dras¬ 
tically changes when the patch size is altered by a factor of 
2. Optimal results are obtained by setting these parameters 
to the values discussed in Section As is shown in Figure 
ED bright filamentary structure is largely insensitive to such 
changes. 

The normalization used in the arctan transform is also 
sensitive to the noise level. When the image is flattened be¬ 
yond the level of the faint filamentary structure, adaptive 
thresholding will return large spurious regions devoid of sig¬ 
nal. We note that such regions can be identified using the 
optimal parameters, however they are eliminated by the re¬ 
moval of small objects (see Figure [^-e). Features near the 
noise level rely on optimal parameter settings to be identi¬ 
fied. The sensitivities discussed above are most noticeable 
in this case and proper identification relies on minimizing 
fragmentation of the region due to noise. 


Normal 

— Convolved 



Figure Bl. Image of the skeleton structures after smoothing 
the Pipe data to the equivalent physical resolution it would have 
at 460 pc. The filamentary structure extracted from the image is 
nearly the same after degrading the resolution by a factor of 2.5. 


APPENDIX B: RESOLUTION EFFECTS 

In this section, we test the derived parameters of the fila¬ 
mentary network with respect to changes in physical reso¬ 
lution. We complete an analysis of the Pipe region, one of 
the nearest targets, at its native physical resolution (0.015 
pc), at the resolution the data would have were the cloud 
located at 460 pc (0.056 pc), and the effect of regridding 
the convolved data back to the original pixel scale. We pro¬ 
cess these versions through the FilFinder algorithm with the 
same parameters adopted throughout this analysis. We dis¬ 
play the identified filamentary networks in Figure |B1| and 
the distributions of the properties in Figure [B^ 

While the network properties and location vary slightly 
with the degradation of the resolution, they are remarkably 
stable despite the poorer resolution. The only noticeable 
change is that there is a tendency for short filaments to be 
merged into longer filaments after image degradation. The 
width of the filaments do not change significantly once the 
data are deconvolved. Overall, the properties do not change 
significantly within the range of resolutions explored, though 
this minimal variation is not expected for data sampling dra¬ 
matically different spatial scales. 

Despite this robustness with respect to resolution, we 
still proceed by degrading all the Flerschel survey data to 


© 2015 RAS, MNRAS 000, 










Filaments 


17 



75 95 99 

Flatten Percentile 


Figure Al. The filament mask and cleaned skeletons from the example image shown in Figure^ are shown with varying patch size 
and the flattening threshold used in the arctan transform. Both parameters have a significant effect on the connectivity of the network. 
Significant features missing in some panels have split into smaller regions, which are rejected by the described criteria for a good filament 
(see The optimal parameter choices are shown in the centre panel. 

a common physical resolution to enable a clear comparison 
between the different systems in the Gould Belt. 
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Figure B2. Comparison of the derived skeleton properties shown 
in |B1| The distributions of the properties are similar under this 
reduction in resolution except that the number of short filaments 
drops by a factor of two. 
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